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ATT ITUni:  DI-TliKMJNATION 


I  ntroduct  ion 


The  orientation  of  the  sensing  axis  of  an  instrument  with  respect  to  a 
fixed  coordinate  system  in  space  is  referred  to  as  its  attitude.  From  the 
attitude  of  a  sensor,  one  may  determine  the  angle  it  makes  with  any  other  vector 
such  as  the  rocket's  velocity  vector,  the  sun  line  vector,  etc.  For  certain 
sensors,  meaningful  interpretation  of  the  instrument's  output  can  be  accomplished 
only  when  the  attitude  of  the  sensor  axis  is  known. 

Various  types  of  gyroscopic  platforms  used  as  attitude  measuring  systems 
on  vehicles  are  discussed.  The  systems  include  the  Attitude  Control  System 
with  two  free  gyros,  the  Miniature  Attitude  Reference  System  with  segmented 
lengths  over  a  5  volt  span  for  the  roll,  pitch  and  yaw  axes,  and  the  Miniature 
Inertial  Digital  Attitude  System  with  a  one-to-one  correspondence  between  angular 
displacements  and  digital  coding  outputs. 

The  next  discussion  centers  upon  processing  procedures  which  produce  final 
attitude  information  from  vehicle  gyro  measurements.  An  analysis  is  shown 
which  relates  any  onboard  probe  vector  P  to  the  loca]  coordinate  system  fixed 
at  the  launch  site.  With  the  orientation  of  P  determined,  the  attack  angle  with 
any  other  attitude  determined  or  predetermined  vector  can  be  supplied.  Often¬ 
times,  the  initial  orientations  of  onboard  sensors  require  corroboration.  A 
study  is  outlined  to  compare  the  phase  relationship  between  onboard  measurements 
and  attitude  predicted  magnetic  field,  lunar,  or  solar  intensities.  This  approach 
can  be  used  for  both  side  and  axial  output,  but  determination  of  phase  angle 
variations  are  restricted  to  well-behaved  areas  of  a  vehicle's  flight. 

To  provide  continuous  final  output  for  the  attitude  system  measurements, 
data  refinement  techniques  arc  introduced  with  the  calculation  of  statistical 
error  values  between  the  predicted  and  measured  output.  Tolerance  levels  are 
set  on  the  measured  input  data  to  ensure  that  all  predicted  data  will  be  within 
a  specified  a  deviation. 


(jyroscojuc  Platforms  Used  as 
A  t' t  i  1 1 1  Jo  Mo  a  s  u  r  i  n  g  Sy  s  t. e m s 

Attitude  Control  Systo in 

Space  Vector  Corporations  Attitude  Control  System  (ACS)  is  comprised  of 
a  Roll  Stabilized  Platform  (RSP)  with  two  free  gyros  which  output  yaw,  pitch 
and  roll  vehicle  motion.  Top  dead  center  (TDC)  of  the  gyro  is  the  reference 
for  all  onboard  probes  (figure  1).  The  roll  and  yaw  are  true  measurements  of 
vehicle  motion  but  since  there  are  only  two  free  gyros,  true  vehicle  pitch  is 
a  function  of  yaw.  further,  if  the  vehicle  yaws  over  85°,  this  will  cause 
loss  of  initial  reference  for  the  system  and  make  attitude  determination 
almost  impossible. 

TDC  Reference 

TDC 


View  from  nose 

Top  Dead  Center  Reference 
figure  1 

Various  models  of  the  ACS  were  used  as  attitude  measuring  systems  on  the 
rockets  analyzed.  All  ACS  models  function  similarly  but  output  differently.  In 
some  cases,  the  pitch  ;uul  yaw  coarse  data  were  calibrated  in  .12  volt  incre¬ 
ments  with  a  linear  conversion  range  of  +90°  corre sponding  to  a  5  volt  span. 

The  pitch  and  yaw  fine  data  over  a  similar  volt  span  had  a  conversion  range  of 
+5°.  The  mil  telemetry  data  had  a  linear  conversion  range  of  +180°  on  a  4.8 
volt  span. 

According  to  figure  2,  pitch  coarse  for  the  ACS  was  decreasing  voltage 
when  pitching  down,  and  a  yaw  right  motion  was  increasing  yaw  coarse  voltage. 

A  clockwise  roll  viewed  aft  to  nose  was  increasing  voltage.  Further,  pitch  fine 
was  increasing  voltage  when  pitching  down,  and  a  yaw  right  motion  was  decreasing 
yaw  fine  voltage. 

i 
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ROLL  ANGLE  INCREASING  CW 


‘LAUNCH  RAIL *0 degrees  ROLL 


ACS  Coordinate  Reference 


Figure  2 


in  other  cases,  only  pitch  coarse,  yaw  coarse  and  roll  gyro  data  were  output. 
The  pitch  data  was  calibrated  over  a  +90°  corresponding  to  a  5  volt  span. 

However,  to  model  the  calibration  information,  discrete  polynomial  expressions 
were  used  over  finite  intervals.  From  0  volts  -  4  volts,  a  linear  model  was 
used  and  from  4  volts  -  5  volts,  a  second  degree  model  best  represented  the 
calibration  information.  The  roll  telemetry  data  had  a  linear  conversion  range 
of  four  i>0°  segments,  each  on  a  S  volt  span. 

Other  situations  existed  whereby  gyro  output  consisted  of  pitch  coarse, 
pitch  fine,  yaw  coarse,  yaw  fine  and  roll  gyro  data.  The  coarse  segment 

was  calibrated  from  -20°  to  +170°  over  a  5  volt  span  with  decreasing  voltage  when 
pi  telling  down.  The  yaw  coarse  segment  was  calibrated  for  a  range  of  +70°  over 
a  3.3  volt  span,  i.c.  ,  .8  to  4.1  volts.  A  yaw  right  motion  was  decreasing  volt¬ 
age  for  the  coarse  yaw  output.  Selected  data  points  were  recorded  for  pitch 
fine  data,  and  the  model  connecting  each  pair  of  points  was  essentially  a  straight 
line.  However,  no  data  calibration  points  were  recorded  for  yaw  fine.  As  pre¬ 
viously  described,  pitch  fine  was  increasing  voltage  when  pitching  down  but  yaw 
fine  was  now  increasing  voltage  when  yawing  to  the  right.  The  roll  segment  was 
calibrated  from  180°  counterclockwise  (CCW)  to  178°  clockwise  (CW)  over  a  4.7 
volt  span  with  voltage  decreasing  for  a  CW  roll  when  viewed  aft  to  nose. 

Another  ACS  gyro  reference  system  consisted  of  pitch  and  yaw  channel  output 
at  +180°  over  a  5  volt  span,  and  roll  channel  output  at  a  360°  segment  over  a  5 
volt  span.  Although  not  used  in  the  data  reduction,  the  pitch  and  yaw  high 
resolution  were  +10°.  The  zero  reference  for  the  gyro  roll  was  the  launch  rail 
(Figure  3).  booking  forward,  yaw  right  i positive,  pitch  up  is  negative  and 
roll  clockwise  is  positive. 


Rail 

ms 


ACS  (ivro  Roll  Zero  Reference 
Figure  3 
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MARS  attitude  system. 


The  Space  Vector  Corporation's  Miniature  Attitude  Reference  System  (MARS) 
provided  gyroscopic  output  for  the  determination  of  rocket  attitude  for  various 
vehicles.  The  roll  axis  coincided  with  the  rocket  axis  at  launch  and  the 
orientation  of  a  gyro  reference  notch  was  used  to  determine  the  location  of  the 
yaw  axis.  The  cross  product  of  the  yaw  axis  with  the  roll  axis  yielded  the 
pitch  axis.  The  axes  and  rotations  of  the  MARS  system  are  identified  in  figure  4. 


Top 


Forward  axis 

MA  RS  Re  f c  r  en c  o  Sy  st  e m 
l;i  guru  4 

The  segmented  lengths  over  a  5  volt  span  for  each  MARS  axis  arc: 

a)  Roll:  90° 

b)  Yaw:  60° 

c)  Pitch:  45°  . 

Segment  identification  of  the  pitch  and  yaw  outputs  for  this  system  required 
continuous  monitoring  of  the  output  to  determine  the  precise  time  of  transition 
from  one  segment  to  the  next.  Hy  keeping,  t  rack  of  the  transitions  and  their 


direction  evidenced  by  .1  j  >  1 1  is  or  minus  I>  volt  ..top  change,  one  could  correctly 
identify  the  segments  presented.  The  roll  output  identification  was  a  5°  +3° 
shorted  segment  in  the  center  of  the  1st  segment  which  appeared  as  a  constant 
output  voltage  for  that  portion  of  the  segment.  This  short  was  repeated  every 
4th  segment . 

MIDAS  attitude  system. 

Space  Vector  Corporation's  Miniature  Inertial  Digital  Attitude  System  (MIDAS) 
provided  gyroscopic  output  used  to  calculate  the  attitude  of  many  rockets.  The 
roll,  pitch  and  yaw  references  were  the  same  as  for  a  MARS  system  but  the  outputs 
differed  considerably.  Whereas  a  MAKS  system  output  in  a  S  volt  span  representing 
some  degree  segment  having  a  predetermined  length,  the  MIDAS  system's  roll,  pitch 
and  yaw  measurements  were  digitized  by  optical  encoders  and  resulted  in  a  one- 
to-one  correspondence  between  angular  displacements  and  digital  codings.  For  the 
digital  coding  variable  n,  the  angular  displacement  4.  in  degrees  was 

*  =  TO 24  '  'H><)  - 

No  bias  values  were  needed  for  conversion  since  the  MIDAS  system  represents 
the  displacement  from  the  uncaged  position  of  the  gyro  which  fixes  the  coordinate 
reference  system.  The  flow  of  attitude  data  for  MIDAS  vehicles  is  simplified 
considerably  when  compared  with  the  flow  of  ACS  and  MARS  vehicles. 


Proctssiai^  Procedures  and  Data  Analysis 


The  oscillogram  containing  vehicle  gyro  measurement s  is  inspected  for  inflight 
calibrations,  data  irregularities,  gyro  malfunctions,  approximate  spin  rates, 
processional  period  and  half  cone  angle.  The  digital  tape  containing,  the  gyro 
data  is  then  unpacked  and  converted  to  degrees  and  plotted.  These  plots  arc- 
then  reviewed  with  results  of  the  oscillogram  inspection  to  identify  obvious  prob¬ 
lems  which  will  adversely  affect  the  data  processing. 

I f  no  gyro  malfunction  occurs  then  the  quick  look  attitude  and  angles  of 
attack  for  the  rocket  axis  are  calculated  based  upon  preliminary  information. 
Should  gyro  malfunctions  occur,  procedures  most  appropriate  to  t ho  available  data 
are  used  to  calculate  the  yaw,  pitch,  and  roll  information  in  the  irregular  area. 

If  the  attitude  measuring  system  is  a  MIDAS,  then  attack  angles  for  requested 
side  probes  are  calculated  at  the  same  time  as  attack  angles  for  the  rocket  axis. 

If  attitude  for  side  probes  is  requested  for  an  ACS  or  MARS  system,  the  roll  data 
is  converted  to  degrees,  checked  for  proper  conversion  and  then  the  preliminary 
side  probe  attitude  is  calculated. 

All  data  irregularities  and  inflight  calibrations  are  filtered  from  the  gyro 
conversions  and  the  results  displayed.  This  output,  if  displaying  a  lack  of 
continuity,  is  curve  fitted  bv  procedures  most  appropriate  to  the  existing  data 
pattern . 

Finally,  if  onboard  probes  measure  magnetic  field,  lunar,  or  solar  intensi¬ 
ties,  it  is  possible  to  corroborate  the  calculated  attitude.  Depending  upon  the 
reliability  of  the  onboard  measurement s ,  a  standard  error  estimate  of  the  attitude 
can  be  included  as  part  of  the  normal  output.  A  flow  of  the  attitude  data  for 
both  ACS  and  MARS  vehicles  is  displayed  in  Figure  S ,  and  a  flow  of  the  data  for 
MI  PAS  vehicles  in  Figure  b. 

Anal ysis. 

Given  the  unit  vectors  X,  Y,  Z.  in  the  directions  of  the  gyro  roll  .  pitch  ( p ) 
and  yaw  (y)  axes,  the  direction  of*  the  longitudinal  axis  of  the  vehicle  can  be 
expressed  as  the  vector  e”: 

e"  =  X  cosy  cosp  +  Y  si  tty  Z  cosy  s  i  np  (11 


Mow  of  Attitude  Data  for  ACS  ami  MAKS  Veiiiclcs 


I  low  of'  Attitude  Data  for  ACS  and  MARS 


Figure  6  (Cont.) 


whore  y  and  p  are  the  true  vehicle  yaw  and  pitch  respectively.  Now  let  cn^, 

c  ,  he  unit  vectors  In  the  directions  of  true  North,  hast  and  the  local  vertical, 
1  c 

fixed  at  the  Launch  site.  Within  this  system,  the  elevation  of  a  vector  is 
determined  as  the  angle  it  makes  with  the  local  horizontal.  The  direction  of 
the  local  vertical  is  positive  elevation.  The  azimuth  of  a  vector  is  the  angle 
between  its  projection  in  the  horizontal  plane  and  measured  positive  clock - 
w L s e  fro m  North  (see  I*  i gu  re  7 ) . 


6  ~  hlevation  of  eM 

r 

6  =  Azimuth  of  e,T 

r 

Local  Last,  North  and  Vertical  System 
Fi gurc  7 


Using  the  direction  coefficients  of  the  gyro  axes  at  launch,  the  unit 
vector-  X,Y,Z  can  be  expressed  as  linear  combinations  of  the  earth-based  system 
t\,  ,  o  ,  er  .  if  we  call  the  coefficient  matrix  B,  then 


We  can  define  a  system  of  orthonormal  vectors  c”,  e”,  e"  at  any  time  in 
flight  as  a  linear  combination  of  the  gyro  axes.  If  A  represents  the  coefficient 
matrix  of  the  gyro  output,  then 


In  particular,  as  described  in  (1) 

=  cosy  cosp 

al2  =  Siny 
a^  =  cosy  sinp 

To  transform  the  time  dependent  o”,  oj',  e"  system  into  the  local  eic> 

er^  system,  form  the  matrix  product  of  A  and  B.  hotting  the  matrix  C  =  AB ,  then 


c.  . 


3 

E  a.  b  .  (i , j  =  l ,2,3) 
.  1m  J 

m=  1  J 


and 


Once  the  attitude  of  the  vehicle's  main  axis  has  been  determined  as  direc¬ 
tion  cosines,  the  elevation  and  azimuth  may  be  derived,  iron  (J),  wo  have 

eM  =  e  cost)  cos<t>  +  e,  cosO  sin<j>  +  c  sin'1 
r  t  S1  r 

c  c  c 


Therefore : 


0  =  sin  1 (c13)  , 


and 


♦  =  tan''  (~)  . 

11 


The  elevation  and  azimuth  of  and  are  similarly  derived. 


In  order  to  extern!  this  approach  to  a  vector  P  lying  in  the  sensing 
direct  ion  of  a  probe  having  any  orientation  on  the  rocket,  P  is  represented 

P  =  o”  cos  A  +  e”  cosp  si  n ;  +  e"  sin^  sin'  , 


where  A  is  the  angle  between  e"  and  P,  and  \i  is  the  angle  between  ey  and  the 
projection  of  P  in  the  plane  of  c’y  and  e y .  furthermore,  if  f  and  are  the 
respective  elevation  and  azimuth  of  P,  it  follows  that 


P  -  c  cost)  cos.'  +  e  cos*  s  in  y  +  e  sin  a 

'!c  P  P  P  P  rc  P 


We  can  also  define  .  -  the  vector  components  of  P  -  by  the  trans¬ 

formation  : 
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where  each  (i ,  j  =  l , 2, 3)  is  a  product  of  c.j  and  the  corresponding  direction 

cosine  of  P.  Using  the  same  technique  as  for  ey,  the  coefficients  of  e<  ,  e  , 

e  in  (a)  and  (4)  are  equated,  readily  yielding  0  and  <}> 
rc  P  P 

finally,  if  V  is  any  vector  associated  with  the  rocket  flight,  the  angle 

between  P  and  V  may  be  derived.  The  unit  vector  in  the  direction  of  V  is  given 

by 


V  = 


+  e  v  +  c  v, 
,•  2  r  3 

c  c 


M 


»  v* 

0  1 
c 


+ 


V  ^ 

r  o 
c 


(5) 


where  the  v.  (i  -1,2,3)  terms  are  the  components  of  V  in  the  directions  of  north, 
cast  and  vertical.  The  resulting  v!  terms  arc  thus  the  direction  cosines  of  Y. 

The  scalar  product  of  (3)  and  (5)  yields 


P*V  =  vj  cos  a  cos<j»  +  v!  cosO  s  i  n  0  +  vl  sinO 

1  P  P  2  p  p  o  P 

Hence,  the  angle  ij;,  between  P  and  V  is  given  by: 

=  cos  ^  f v !  cosO  cosO)  +  vi  cosO  sin<i>  +  v’  sinO  ) 
L  p  p  2  P  P  3  )V 
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If  the  unit  vector  V  is  replaced  by  an  offboard  unit  vector  representing 
magnetic  field,  lunar,  or  solar  Intensities,  then  this  offboard  vector  can 
take  the  form 

H  =  e.  h-  +  e .  h0  +  e  h„ 

0  1  f  2  T  ^ 

c  c  c 

The  geocentric  coordinates  together  with  the  launch  data  determine  the  direction 
cosines  h^,  h9,  h^  of  the  offboard  unit  vector  li.  The  pitch  angle  between  11 
and  the  onboard  sensor  P  in  the  direction  of  the  magnetometer,  solar  sensor, 
or  lunar  sensor  is 

6^  =  cos  *(H*P) 

However,  the  pitch  angle  can  also  b^.  computed  directly  from  onboard 

measurements  as  the  cos  of  the  ratio  of  the  intensity  in  the  P  direction 
to  the  total  intensity.  The  pitch  angles  and  ft,  are  then  compared  to  confirm 
the  attitude  data,  and  check  for  pha  *  angie  variations. 

Data  refinement. 

Due  to  irregularities  oft^n  present  in  the  recorded  data,  filtering  pro¬ 
cedures  are  modified  or  developed  as  required  to  provide  continuous  final 
output  for  the  attitude  determination  system.  The  primary  routines  used  are 
Fourier  Series,  nth  degree  polynomial  approximations,  and  moving  average 
methodology. 

To  provide  smooth,  continuous  and  acceptable  attitude  information  for 
vehicles  when  they  have  established  angular  preccssional  velocity  and  a  half 
cone  angle,  it  is  possible  to  generate  gyro  yaw  and  pitch  information  by  means 
of  the  Fourier  Series  expansion.  During  this  well-behaved  area  of  a  particular 
vehicular  flight,  the  coarse  and/or  fine  yaw  and  pitch  data  can  be  predicted 
by  the  Fourier  Series 

a  00 

o  ,  mrx  ,  nnx 

y  =  —  +  >.  (a  cos  -7—  +  b  sm  ) 

7  2  .  n  L  11  L 

n=  1 


Tins  procedure  has  proven  successful  due  to  the  periodic  nature  of  the  data. 
However,  for  most  rocket  flights,  quick  convergence  of  the  Fourier  Series  is 
obtained  with  an  approximation  i  nco  rporat  ing  a  linear  term 

V  ~  f.  +  f,t  +  SinMt  +  f.  COSut 

1  2  o  4 

A  preliminary  angular  velocity  u-  is  selected  from  a  study  of  the  oscillograms 
of  the  raw  data,  and  is  further  refined  by  an  option  within  the  fitting  routine. 

During  regions  in  which  a  vehicle  is  not  coning,  yaw  and  pitch  outputs  are 

separated  into  discrete  time  intervals.  These  intervals  are  then  curve  fitted 

t  h 

with  polynomials  up  to  the  9th  degree  when  necessary.  The  n  degree  polynomial 
routine  calculates  an  RMS  value  between  measured  and  predicted  data  for  each 
of  the  polynomial  approximations.  The  minimum  RMS  value  determines  the  degree 
of  the  polynomial  approximation  to  be  used  in  the  specified  interval. 

Moving  average  methodology  has  also  been  used  to  correct  for  random  noise. 
Variations  have  included  n  by  m  weighted  procedures  to  follow  true  data  patterns. 
This  approach  is  used  when  continuous  but  irregular  measurements  are  available, 
and  smoothing  is  required. 

Regardless  of  the  refinement  technique  employed,  a  standard  error  estimate 
1  is  computed  between  measured  and  calculated  values  to  compute  an  acceptable 
confidence  interval.  This  confidence  interval  is  used  to  select  the  optimum 
modeling  approach  among  those  under  consideration. 

V ehicle s  pro c essed . 

Using  the  procedures  described  in  this  report,  20  attitude  reports  for 
various  vehicles  have  been  prepared  during  the  reporting  period  from  1  January, 
1979  to  31  March,  1981.  A  listing  of  completed  vehicle  attitude  appears  in 
fable  1.  Final  outputs  consists  of  launcher  and  vehicle  referenced  plots  and 
listings  of  selected  parameters.  A  final  attitude  tape  containing  all  available 
time  referenced  outputs  is  also  generated. 


1( 


TRAJECTORY 


t 


lntroduct ion 

Raw  radar  or  telemetry  tracker  recorded  polar  coordinate  data  for  a  rocket  is 
normally  available  on  magnetic  tape.  However,  such  tapes  are  not  initially  compat¬ 
ible  with  the  Control  Data  Corporation  (CDC)  60  hit  word  format.  To  enable  process¬ 
ing  of  the  data,  the  tape  must  be  unpacked  and  rewritten  in  a  CDC  compatible 
form.  Once  unpacked,  the  data  is  examined  for  excessive  noise  and  prolonged 
missing  data  intervals.  If  either  condition  exists,  appropriate  data  plugging 
procedures  are  initiated  prior  to  a  geocentric  coordinate  conversion  of  the  radar 
or  tracker  referenced  slant  range,  azimuth,  and  elevation  data. 

In  the  event  multi  radar  and/or  tracker  information  is  available,  conver¬ 
sions  of  the  available  data  to  any  desired  reference  location  is  possible.  This 
transformation  enhances  the  reliability  of  the  filtering  procedures  when  exces¬ 
sive  noise  or  missing  data  in  selected  intervals  cannot  be  modeled.  in  such 
cases,  reliable  positional  data  from  appropriate  radars  and  trackers  are  merged 
to  form  a  continuous  data  base. 

When  an  adequate  model  for  the  vehicle  is  available,  normal  data  filtering 
procedures  consist  of  optimizing  the  launch  conditions  to  the  start  of  the 
filter  limits  by  back  filtering  the  data  to  launch.  When  a  model  is  not  avail¬ 
able,  then  a  double  filter  is  normally  performed.  This  procedure  back  filters 
the  positional  data  to  only  the  start  of  the  filter  limits.  in  both  cases,  either 
integration  or  data  reproducing  options  are  available  beyond  the  filter  limits. 

Final  outputs  for  the  rocket  trajectory  system  consist  of  filtered  position¬ 
al  values  referenced  to  the  launcher  and  the  radar.  Residual  data  during  the 
back  filter  period  is  also  available. 

Radar  and  Tracker  Systems 

Several  tracking  systems  have  been  employed  during  this  reporting  period. 

A  brief  description  of  the  primary  unpacking  routines  for  these  systems  follows: 

1.  RTMIN  -  used  for  real  time  radar  data  from  Whi  te  Sands.  Positional 
values  are  launcher  referenced. 

2.  REFORMS  -  used  to  reformat  radar  referenced  data  received  from  the 
Wallops  Might  Center  (WFCJ. 
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3.  REFORMN  -  used  to  reformat  PNA  tracker  data  from  Poker  Flat  Research 
Ran ge . 

4.  REFORMO  -  used  to  unpack  telemetry  tracker  data  from  the  Oklahoma 
State  University  (OS U j  Tradat  system. 

5.  3x80  -  used  to  reformat  single  radar  tracking  data  from  White  Sands. 

6.  3xK0  MULT t  -  used  to  reformat  multi  radar  tracking  data  from  White  Sands. 

Additional  features  were  incorporated  into  the  3x80  multi  unpacking  routine 
in  order  to  more  adequately  compare  the  vehicle  trajectory  results  from  each  radar 
employed.  Residuals  between  x,y,z  and  x,y,z  are  calculated  and  plotted  for  desired 
pairs  of  radars.  These  displays  are  used  as  aids  in  judging  which  radar  is  most 
rcl i able . 

The  unpacked  data  is  written  in  16x6  blocks  with  radar  referenced  time, 
slant  range,  azimuth,  and  elevation  data  being  the  normal  configuration  in  each 
16  word  record.  In  the  event  excessive  noise  or  missing  data  exists  for  any  of 
the  parameters  in  this  file,  data  plugging  and  modeling  procedures  are  used  to 
rewrite  the  16x6  data  file  in  a  form  compatible  with  the  remainder  of  the  traject¬ 
ory  processing  system. 


Processing  Procedures 

Should  excessive  missing  data  or  noise  exist  in  the  unpacked  data,  depend¬ 
ence  on  a  single  modeling  procedure  for  the  ill  behaved  areas  are  often  insuffi¬ 
cient  to  produce  reliable  results.  Hence,  prior  to  geocentric  conversion  of  the 
vehicle  positional  data,  care  is  taken  to  provide  a  most  complete  and  reliable 
estimate  of  the  slant  range,  azimuth,  and  elevation  information  over  time. 

Pat  a  nodi ficat i on . 

To  correct  extremely  random  data  patterns  requires  either  establishing  a 
confidence  band  for  acceptable  values  or  generating  a  time  series  model  through 
the  data.  In  the  case  of  establishing  a  confidence  band,  a  weighted  moving 
average  of  the  data  in  selected  intervals  is  created.  The  residuals  between  the 
raw  and  moving  average  data  are  then  calculated  and  a  standard  error  estimate 
a  is  computed.  A  P.r>!.  confidence  interval  for  acceptable  measurements  is  normally 
established.  Any  time  related  raw  positional  value  outside  this  interval  is 
deleted  from  data  base  of  acceptable  values,  and  is  replaced  by  the  weighted 
moving  average  value. 
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When  existing  data  patterns  are  readily  obvious,  but  stiJl  randomness 
exists  over  extended  periods,  polynomial  regression  procedures  are  an  alternate 
data  modification  method  to  the  establishment  of  confidence  bands  for  acceptance 
or  rejection  of  the  data.  In  such  cases,  higher  ordered  polynomials  are  succes¬ 
sively  passed  through  the  data  being  examined.  Residual  error  between  the  raw 
and  calculated  data  is  determined  by  the  RMS  value.  The  minimum  RMS  is  the 
criterion  for  selection  of  the  appropriate  polynomial  model  . 

To  fill  missing  data  intervals  (greater  than  10  seconds),  poly¬ 
nomial  modeling  is  used  only  when  data  patterns  are  quite  apparent.  In  the 
event  of  irregular  motion  when  multi  radar  information  is  available,  reliable 
positional  data  from  appropriate  radars  and  trackers  are  merged  to  form  a  con¬ 
tinuous  data  base  prior  to  the  data  modification  methodology.  This  procedure 
requires  a  coordinate  transformation  of  the  radar  or  tracker  information  to  the 
desired  reference  point. 

Transformation  of  positional  informat  ion. 

Given  two  or  more  sets  of  positional  information  for  a  given  vehicle  from 
different  reference  points,  merging  of  selected  intervals  from  each  data  base 
may  be  accomplished  in  order  to  generate  a  more  continuous  and  reliable  file  of 
positional  information.  Assuming  the  geodetic  latitude,  longitude,  and  height 
above  the  Earth!s  surface  of  the  given  reference  point  to  be 

Oj  :  latitude 
:  longitude 
:  height, 

and  the  reference  of  the  desired  system  to  be  U?,  and  II  v  the  goal  is  to 

transform  the  given  slant  range,  azimuth,  and  elevation  data  to  the  new  location 
The  raw  positional  data  shall  be  designated  as 

SR^  :  slant  range 
Aj  :  azimuth 
H1  :  elevation  , 

and  the  desired  information  as  SR0,  A,,  b^. 


[n  the  geocentric  system,  the  vector  to  the  target  vehicle  is  defined  as 
R  =  r  ( i  cos-’^cos^j  *■  jcos  ^sin^  +  ksint^)  (6j 


where 


r  =  \/x“  +  y  + 


lz+roJ‘ 


x  =  SR^  cosily  sinAj 


y  =  SR^  cosily  eosAj 


2  =  SR1  sinf^ 


r  =  637cS.  185  -  21.30  sin"o  . 
0  £ 


'Hie  angles  n  and  ;»  are  the  geocentric  latitude  and  longitude  of  the 
target  vehicle.  The  geodetic  latitude  is  computed  as 

tanh  =  l.OOhfc 74  tana  . 

£  h 

Furthermore,  given  t he  geocentric  latitude  and  longitude  ^  of  the  initial 
radar/ tracker. 


sineL  = 


y  cosO.  f  (z  +  r  J  sin  ' 

1  o'  1 


-y  sine  sinL  +  x  cos<L  +  (z  +  r  )  cosh  sin<L 

tan  r  =  - - - - - t - - - - - i - -  . 

L  -y  smO  cos-;*.  -  x  siniK  +  l  z+r  )  cosO,  cost), 

11  loll 

From  (<>)  ,  the  components  of  the  geocentric  referenced  position  update 


vector  are 


XV  =  r  cos  ^cos+j 


YV  =  r  cosO  sin*}). 


Z\‘  =  r  sinO. 


The  problem  is  now  to  find  a  transformation  matrix  k  to  transform  the 
original  radar/ tracker  site  data  to  the  desired  site.  To  calculate  SR^  ,  A0, 


and  li^  for  the  offset  site  where 


AN2  *  SR^  cosE^  cosA,  (<jj 

AV2  »  SR2  sinE9  , 

we  must  first  express  the  local  system  components  AE2,  AN2 ,  AV2  at  the  offset 
site  in  terms  of  available  information.  In  matrix  form,  let 


Computation  for  the  Euler  matrix  K  and  offset  vector  X  ,Y  ,Z  follow. 

1  os  os  os 

In  the  i,j,k  geocentric  system,  let  P^  be  the  position  of  the  initial  rad 

tracker  and  P2  the  position  of  the  desired  offset  site. 

Pi  =  (R^+Hj)  (icoso^  cos*  +  jeoso  sin<f^  +  ksina  J 

P?  =  (icos^2  cos^?  +  jeoso,  sin«^  +  ksim  ,) 

The  radii  R..  arc  computed  as 

R.  =  6378.185  -  21.39  sin2o  ,  i=l,2 

1  *i 

where  0^  is  the  geodetic  latitude  of  the  respective  radar/ t racker. 

The  offsets  from  the  geocentric  to  the  offset  location  on  the  Earth* 
surface  are 

Xos  =  (R2+I!2)  cos')2  cos4-u 

Yos  =  (R2+1I2J  cos "2  sin<!;|) 

Zos  -  («2.ll,)  ,lm.2 

where  is  the  longitude  offset  from  the  initial  radar/tracker  to  the  offset 
site. 

The  distance’  vector  I)  (Tiyuri  8)  between  Pj  am!  P(  i  ik  lin-  J  a  - 


1)  -  i[(K^+H^)  cos  >  cos:,  (Rj+Hj)  cosf'^  coss^) 
+  j[(R„+i!,)  cos  ,  sin^  -  (R1+H1)  cosu^  sino^] 
+  k[(R7+il0)  s  i  n  *  * .,  -  (1^+lljJ  sin  }1  - 


Vector  Difference  Between  Rauars 
Figure  8 

Simplifying  equation  (10)  in  terms  of  direction  cosines, 

D  =  - - -  (a  i  +  a  i  +  a  k) 

iv,  i  1  2  3 

where 

|P2-pi!  =  |/af  +  a 2  ~  aj  . 

In  the  i,j,k  system,  the  distance  vector  can  so  be  expressed 
\)  ~  icosO^  cos<}>^  +  j co s 0 p  sin<^  +  ksinU^  . 

Equating  the  two  expressions  for  p. 


i  1 


(10j 


These  rotation  angles  0  and  <|ip  relate  the  local  componen t  s  at  the  do-;  i  red 
offset  site  to  the  original  site  information  as 


-sineD  cos*n 

^ COS Op  COS^p 


-sin Op  sin*D 
cosi)p  sin<J»p 


VV  -  Y 


Equating  systems  (9)  and  (11),  the  desired  parameters  at  the  new  site  are 

SR2  =  \/aE22  +  AN22  +  AV2~2 
A2  =  tan-1  (AE2/AN2) 

E2  =  sin'1  (AV2/SR2) . 

Corrections  for  extended  missing  data  intervals,  extreme  data  irregularity, 
and  merging  of  multi  radar  information  having  been  made,  the  trajectory  data 
flow  now  continues  with  the  geocentric  conversion  of  positional  values  [see 
Figure  y) , 

Using  corrected  radar  or  tracker  slant  range,  azimuth  and  elevation  data  the 
next  phase  of  trajectory  determination  calculates  geocentric  position  and 
velocity  information  as  well  as  geodetic  altitude  of  the  vehicle.  Using  a  four 
point  cubic  polynomial,  further  adjustments  to  the  raw  positional  data  can  be 
made  if  necessary.  Listings  and  displays  arc  reviewed  for  data  stability  and 
completeness  before  continuing  with  the  filtering  procedures  to  produce  final 
trajectory  data. 

Fi  ltering. 

Depending  whether  or  not  a  model  for  vehicle  dynamics  exists,  the  corrected 
data  base  consisting  of  time,  slant  range,  azimuth,  and  elevation  data  is  passed 
through  an  appropriate  filtering  procedure.  When  a  model  does  exist,  a  launch 
optimi  zing  procedure  is  normally  performed.  In  this  case,  filter  limits  arc 
set,  and  the  launch  conditions  are  optimized  to  the  start  of  the  filter.  When  a 
model  is  not  available,  a  double  filter  which  back  filters  the  positional  data 
to  only  the  start  of  the  filter  limits  is  the  normal  procedure.  Filter  limits 
are  normally  selected  as  the  interval  of  well-behaved  performance,  and  in  all 
cases  include  data  above  50  kms. 


Regardless  of  the  filtering  procedure,  both  integration  and  data  reproducing 
options  arc  considered  to  optimally  complete  the  vehicle  trajectory  at  low  alti¬ 
tudes  (below  40  kmsj.  Often,  increasing  adjustments  to  the  base  angular  variance 
are  necessary  for  relatively  smooth  launcher  referenced  velocities.  Examination 
of  the  plotted  results  and  analysis  of  the  data  structure  normally  suggest  the 
required  adjustments  at  this  phase. 


Vehicles  processed . 

Using  the  procedures  described  in  this  section  of  the  report,  18  trajectory 
reports  have  been  prepared  during  the  reporting  period  from  1  January,  1979  to 
31  March,  1981.  A  listing  of  completed  trajectories  appear  in  Table  1,  Final 
outputs  consist  of  launcher  and  radar  referenced  plots  and  listings  of  selected 
positional  parameters.  A  final  trajectory  tape  containing  all  available  time 
referenced  filtered  outputs  is  also  generated. 


Tabic  1 

COMPUTED  ATTITUDE  AND  TRAJECTORY  RUNOUTS 


Vehicle 


At  t  i  tud  v  Tr  a  j  e  c  t  o  r  y 


A2  4. 609-2  (MSMP) 

Booster  x  x 
Target  x  x 
Sensor  x  x 


IC819.08-1 

Nose  cone  x  x 

Payload  x  x 


IC807. 15-1 
IR807.S7-1 
A10. 802-1 
A10. 802-2 
AO  7. 712-1 
A08. 705-2 
AO 8. 708-1 
A12.9A1 
A12.9A2 
A24.7S2-1 


A18.805 

MIM  x 

PIM  x 

A31 . 702  x 

A04.703  x 

PLUMEX1  x 

PLUMEX2  x 

WS834.28-2  x 

WS834.28-3  x 

A10. 901-1  x 

A51.970  (EXCEDE  II)  x 


x 


x 


X 


* 


Due  to  special  data  problems  encountered,  efforts  are  continuing 
in  the  determination  of  reliable  positional  information  for 


HXCHDH  I  I . 
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ANALYSIS  TliUINlQUFS  APIMJ  li!)  TO  BALLOON- BORN  I-  MOSAIC 
I  NTERF’LROMITLR  AND  RADIOMETER  MEASUREMENTS 


Computational  techniques  developed  for  the  reduction  and  analysis  of 
interferometer  and  radiometer  data  from  the  Balloon  Altitude  Mosaic  Measure¬ 
ments  IBAMM)  program  arc  presented.  limphasis  is  given  to  the  reduction  and 
analysis  techniques  developed  for  mosaic  detector  interferometer  measurements. 
The  presentation  includes  discussion  of  techniques  to  quality  check,  edit  and 
simulate  radiance  values.  Statistical  analyses  performed  on  simulated  radiance 
data  is  described  and  examples  of  raw  and  edited  data  and  statistical  param¬ 
eters  are  given.  A  brief  discussion  of  radiometer  reduction  and  analysis 
techniques  is  presented  which  includes  data  editing;  filtering  and  decimation; 
statistical  parameter  determination;  autocovariance  estimation;  and  Fourier 
Rower  Spectral  estimation  of  radiance. 

Int  roduct ion 

This  djsctission  concerns  the  digital  signal  processing  techniques  developed 
to  analyze  and  reduce  radiometer  and  interferometer  measurements  into  radiance. 
ITte  use  of  interferometer  measurements  allowed  the  determination  of  the  radiance 
values  over  any  wavelength  band(s).  These  radiance  profiles  are  then  statisti¬ 
cally  analyzed.  The  initial  major  concerns  in  developing  the  analysis  were  the 
magnitude  of  the  error  that  would  be  added  to  the  data  base,  because  of  the 
mathematical  limitations  of  the  techniques  utilized  and  computer  time.  These 
concerns  were  minimized  by  the  fact  that  the  noise  inherent  to  the  experiments 
dominated  those  that  would  be  introduced  by  the  mathematical  models.  The  pro¬ 
cedures  and  computer  coding  were  developed  to  minimize  the  amount  of  computer 
time  and  place  some  limitation  on  the  computer  storage  used. 

Since  the  analysis  required  the  calculation  of  many  discrete  Fourier 
Transforms,  it  was  essential  to  develop  a  fast  version  of  the  Fast  Fourier 
Transform.  (The  last  Fourier  Transform  is  an  algorithm  for  rapidly  calculating 
Discrete  Fourier  I  runs  forms .  )  The  routine  developed  minimizes  the  time  spent 
deciding  which  "butterfly"  (butterfly  refers  to  a  logical  procedure  used  mar/ 
times  in  the  Fast  Fourier  Trans  lorm  technique)  to  calculate  next.  Other 
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techniques  such  as  suggested  by  Morris1,  are  available  but  the  one  used  was 
judged  the  fastest  given  the  computer  storage  limitations. 

Data  characteristics. 

The  sample  rate  for  the  radiometer  data  base  (one  sample  every  .0073 
seconds)  was  such  that  very  little  aliasing  (ambiguity  in  deciding  at  which 
frequencies  the  Fourier  Power  exists)  occurred.  Further,  there  was  negligible 
Fourier  Power  near  the  Nyquist  Frequency  (the  highest  frequency  inherent  in 
the  Discrete  Fourier  Transforms). 

The  interferometer  measurements  consisted  of  a  ’’scan”  every  At  seconds. 
Sixteen  simultaneous  interferograms  were  generated  for  each  scan  (only  the 
central  4  detectors  were  used  in  a  rapid  scan  mode).  Rach  detector,  however, 
is  processed  individually.  The  focal  plane  consisted  of  a  4^4  mosaic  array  of 
detectors,  hence  for  the  two  scan  rates  At,  was  either  .25  seconds  or  .0575 
seconds.  A  scan  consisted  of  measurements  of  the  i nterferogranis  starting  at 
maximum  optical  retardation  and  moving  to  the  other  extremum.  The  next  scan 
repeated  the  process  heading  in  the  opposite  direction.  Limitations  in  the 
equipment  meant  that  the  time  interval  between  every  other  scan  was  constant. 

The  time  interval  between  consecutive  scans,  however,  was  a  function  of  the 
direction.  Similarly,  the  mechanical  position  corresponding  to  the  zero  optical 
retardation  for  every  other  scan  was  reproducible  but  differed  for  the  two  scan 
di rections. 

The  sample  rate  for  measurements  for  each  scan  was  high  enough  to  prevent 
aliasing.  At  was  short  enough  that  changes  in  the  interferometer  pattern 
during  a  scan  could  be  ignored  and  further,  that  only  a  small  amount  of  alias¬ 
ing  occurred  when  analyzing  a  set  of  scans.  When  events  occurred,  however.  At 
was  not  short  enough  to  assume  that  there  was  negligible  Fourier  Power  near 
the  Nyquist  Frequency  of  a  set  of  interferometer  scans. 

Filtering  and  decimation. 

The  processing  of  the  edited  radiometer  data  takes  large  annum ts  of 
computer  time.  In  order  to  speed  the  processing,  the  highest  frequency  of 
interest  was  determined.  Given  this  highest  frequency  (F^)  a  low- pass  filter 


was  designed  and  applied  to  the  data.  The  filter  used  was  designed  to  fit  the 
following  frequency  response  curve  in  a  least  square  sense. 
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Nyquist  Frequency 
Figure  10 


Here  F  was  arbitrarily  chosen  approximately  equal  to  1.1  1  .  This  filter  was 

a  symmetric  finite  impulse  response  filter.  hnough  coefficients  were  used  so 
that  the  response  curve  nearly  coincided  with  the  above  model. 


This  filter  was  applied  in  frequency  space  using  the  convolution  theorem. 
The  filtered  data  base  was  then  decimated  (every  Nth  point  used)  so  that  the 
new  Nyquist  frequency  of  the  new  data  base  was  F^,  the  start  of  the  stop  band 
of  the  discrete  filter  (Fourier  Power  above  F^  approximately  0.).  This  pro¬ 
cedure  prevents  any  significant  aliasing  of  the  decimated  data  base.  The  amount 
of  decimation  was  generally  not  a  whole  number  of  data  points.  In  order  to  find 
the  value  between  filtered  data  points,  as  specified  by  the  decimation,  linear 
interpolation  was  used.  The  radiometer  values  after  filtering  did  not  have  a 
large  amount  of  Fourier  Power  at  frequencies  near  F^.  This  meant  that  linear 
interpolation  added  negligible  error  to  the  analysis. 


Data  Processing  Approach 

The  original  radiometer  data  base  occurred  at  approximately  equally  spaced 

time  intervals.  (bus  we  had  data  pairs  R.  ,  t^  where  R.  is  the  radiance  reading 

at  time  t..  In  the  original  data  base  some  of  the  R.  values  were  unavailable 
)  ’  } 

.aid  others  appeared  to  be  unrealistic.  ihrec  reasons  for  unrealistic  ('bad”) 
measurements  arc  jxjssible.  One  is  that  the  data  is  ’’stationary”  and  the 


measurement  is  real.  The  odds  of  this  happening  are  negligible  (less  than 
one  in  ten  thousand)  with  this  analysis.  A  second  explanation  is  that  the  ’'bad” 
points  are  "gliches"  in  the  measurements  due  to  problems  in  the  collection, 
transmission,  digital  conversion  or  processing  of  the  data.  The  procedure 
attempts  to  remove  these  "gliches".  The  third  explanation  is  that  the  ’’bad" 
points  are  actual  "events".  in  order  to  see  If  "events"  are  removed  by  the 
editing  procedure  the  effect  of  editing  on  measurements  of  known  "events"  was 
checked.  The  method  selected  for  editing,  the  data  involved  first  differences 
and  as  proceeded  follows: 

1)  Calculate 

A  »  EUR.-Rp/U-i) 

where  the  sum  is  for  all  pairs  of  points  satisfying  the  following 

criteria.  R.  and  R.  are  available,  j>i  and  i-i  the  minimum  such 
j  l  * 

value.  N  is  the  number  of  terms  in  the  sum. 
s 

2)  Similarly  calculate 

MS  =  (((R.-Kp/U-i))2) 

S  J 

3)  Remove  "bad"  R^  values  until  al l  the  remaining  points  satisfy  the 

criteria 

|  ((R.-R^/U-i))  -  A |  VMS*S 

where  S  is  a  parameter  set  to  4. 

For  the  radiometer  this  editing  had  no  discernible  effect  on  known  "events" 

After  "bad"  points  were  removed  from  the  radiometer  data  base  the  missing 
radiance  values  are  filled  in  using  linear  interpolation.  As  long  as  there 
are  no  long  time  intervals  where  the  data  was  missing  and  there  is  negligible 
Fourier  Power  near  the  Xyquist  frequency  linear  interpolation  was  sufficient. 

Radiometer  or  Radiance  Calculated  from  Interferometer  Intensity 


Statistical  estimates  were  used  as  tools  to  validate  or  formulate  theories 
of  how  the  actual  data  behaved.  The  estimates  used  in  our  processing  are 
exceedance  (FIX),  Probability,  Autocovariancc  (All),  and  a  Power  Spectral  Density 


E  s tj w t  i  excco da n c e . 

Exceedance  is  defined  as  EX(X)  =  l-E(X)  with  i  ( X)  the  cumulative  probability 
fun  ct ion . 

The  following  method  of  estimating  the  exceedance  of  (RAD)  (radiance  values j 
was  used: 


1)  Estimate  the  mean  as  AVI:  =  ~  >:  V. 

N  L 

where  V.  are  the  K;M>  values  and  N  is  the  number  of  measurements 

i 

2)  Estimate  the  standard  deviation  as 
<  --  y>:(Y.  -.\vi.)"/iN-l  > 


xi  Divide  the  interval  x  +  AVI:  and  3.t  +  AVI:  into  400  equally  wide  bins. 
Add  one  additional  bin  for  all  values  less  than  -3*  +  AVE.  Order 
these  bins  according  to  their  radiance  values.  Place  and  count  the 
V.  values  into  their  appropriate  bin  (BL\\)  where  j  =  l  ,2 , .  . .  ,401) 

4 )  The  exceedance  at  the  high  edge  of  the  jth  bin  E(RAlh)  is  then  esti¬ 
mated  as 


r.X(KAi).)  i  - 

1 '  iN 


Es  t  i  ina t  i  nj£  jrmb  ab  i  l  i  t  y . 

The  probability  of  radiance  occurring  in  a  particular  interval  was  esti¬ 
mated  using  the  first  three  steps  in  the  exceedance  estimation,  with  a  smaller 
number  of  bins.  These  values  could  also  have  been  calculated  from  estimates 
of  the  probability  density  function,  which  could  be  estimated  by  several 
di fferent  methods 


Es t  imaf  ing  autocovar/ancc. 

The  autocovariancc  function  of  a  set  of  discrete  data  0\)  is 
A I  i  (  i)  -  E((Y  -  V  I  ( V  .  -V  )  i 

v  i  i  + 1 


whore  ^  is  the  mean  of  Y.  and  I  is  the  expected  value  and  ;  is  known  as  the  lag. 


Two  methods  of  estimating  the  autocovn r iance  of  the  radiance  were  used. 

The  estimates  are : 

Estimate  (1) 

AU(1)(t)  =  JjZ((Y.-AVF.)  ( Y.  +  (-AVfi) ) 
which  is  a  biased  estimate  of  the  autocovariance  and  estimate  ( 2 ) 

au(2)(t)  =  T(Yj'AVii)  (VcAVR)) 

which  is  an  unbiased  estimate  of  the  autocovariance.  It  should  be  noted  that 
AVK  could  be  replaced  by  some  higher  order  curve  fit  of  the  data  and  that  other 
reasonable  estimates  of  AU  are  possible. 

listimat  ing  Pow cr  Spectral  l )en s i  t y 

The  integral  of  the  Power  Spectral  Density  between  some  frequencies  fj  and 
f?  is  the  energy  found  in  that  frequency  interval.  The  total  energy  is  assumed 
to  equal  E((Y^-Y)^).  Where  I:  is  the  expected  value  and  Y  -  E(Y.). 

Three  methods  of  estimating  the  Power  Spectral  Density  of  RAD  values  were 
used.  In  all  these  methods  the  final  results  were  normalised, 
fs  1 

Let  P ^  be  the  unnormn  1  i zed  Power  Spectral  Density  estimates  ot  RAD  where 

P^s^  is  the  s ^  type  estimate  of  the  Fourier  Power  in  the  frequency  If)  interval 
( i - 1 ) *A  f<f  -  —•  •  i  *A  f .  I  le  re  A  f  is  t  he  reso  1  ut  i  on  f  requeue  y  o  f  t  he  Four  i  e  r  Pow  c  r 
estimates.  Then  multiply  P.  ^  by  a  constant  such  that 

,  ,  L  N 

Nflsj  Af  >:  P.  J  i  }  (Y.-AVli)" 

.  .  l  N  .  i 

i-l  J=1 

where  L*Af  is  the  Nyquist  Frequency.  The  estimates  P.^  =  P.^  are  the 

normalized  Fourier  Power  Fist  i mates. 

This  normalization  was  used  rather  than  attempting  to  estimate  the  weight¬ 
ing  factor  of  any  apod izat ion  (windowing)  or  smoothing  performed  to  find  the 
Fourier  Power  Estimates. 


The  three  methods  used  to  estimate  the  unnorma  1  i  ::ed  Fourier  Power  were: 


Met  ik>d  i  -  Uaw  Periodogram  Pst  imat  es 

The  per  iodog  nun  estimate  involves  finding  discrete  Fourier  Coefficients 
of  lv  \  i )  and  using  the  square  of  these  coefficients  as  the  unnormal  i  ted  Fourier 
Powt r  estimate. 

Method  II  --  Lstimate  from  Bias  Autoeo variance 

This  method  involved  truncating  to  form  B(kj,  B(i)  is  then  apoui:: 

using  a  triangular  window  (Bartlett  Window),  although  other  apod i tat  ions  could 
have  been  used.  flic  unnormal  iced  Fourier  Power  estimates  were  the  cosine  sene 
coefficients  of  the  resultant  values. 

Method  111  -  lstmiate  from  Unbiased  Autocovariance 

(2) 

This  method  is  exactly  the  same  as  Method  II  but  applied  to  AIJ  (;). 

In  general  the  three  methods  usually  lead  to  approximately  equal  values. 

It  should  be  noted  that,  there  are  other  methods  of  estimating  the  Fourur 
Power  by  fitting  the  data  to  models.  Some  of  the  more  popular  such  methods 
arc  known  as  Maximum  iaitropy.  Maximum  Likelihood  and  the  Autoregressive  Moving 
Average  methods  of  Fourier  Power  P.st  imat  ion. 

'I he  estimates  calculated  by  these  three  methods  might  be  improved  by  pre- 
win  ten  ing  (filtering  the  data  to  better  approximate  white  noise)  the  data  and 
then  post-darkening  (correcting  the  Power  Spectral  Fsti mates  for  the  transfer 
function  of  the  pre -whiten ing  filter)  the  results. 

Pd i ( i ng  I nterfe rome ter  Data 

from  each  scan  of  mt  eriVromet  i-r  data  one  value  of  simulated  radiance  for 
a  particular  f requeue)  baiuu  ■■)  was  determined.  I  he  sample  rate  of  the  scans 
is  such  that  it  was  adv  i  • 1 1 » 1  *•  remove  as  few  scans  as  feasible. 

An  interfemgram  oft.n  tor-*  t ; .  c  f^cn  that  tin  ends  of  the  data  occurred 
at  approximate  nulls.  in  i  •  *  *  t  u  *  «  i «  me  noils,  the  data  tended  to  oscillate  with 
larger  and  large]'  amplitude  toward  i  in  coaler  of  the  in t erferogram.  l\iier.  the 
data  had  nulls  at  its  onb>,  a  special  editing  technique  was  used  to  detect  any 
unrealistic  i "had")  points  i  n  tn<  i«  of  the  nulls.  This  technique  accepted 

measurements  which  either  fell  within  ..  band  around  t  lie  \alue  of  the 


approximate  null  or  had  a  neighboring  point  with  approximate1])  the  same  value. 
Whenever  a  "bail"  point,  was  detected  it  was  set  to  a  neighboring  value. 


The  measurements  excluding  any  nulls  was  then  arbitrarily  divided  into 
five  equal  segments.  hach  of  those  segments  was  edited  sepa rat  cl v .  The  edit¬ 
ing  of  one  of  these  segments  involved  finding  the  data  point  in  this  segment 
which  caused  the  largest  magnitude  of  first  difference  of  the  data.  This  ;>oint 
was  deemed  "bad"  if  the  magnitude  of  the  first  differences  associated  with  it 
was  more  than  twice  the  magnitude  of  the  neighboring  first  differences.  When¬ 
ever  a  point  was  "bad",  it  was  replaced  using  linear  interpolation.  After 
replacing  a  "bad"  point,  this  algorithm  was  repeated  for  the  same  data  segment. 
When  the  point  checked  was  not  "bad",  the  procedure  was  begun  on  the  next  data 
segment  (if  any).  If  there  were  an  unacceptable  number  of  "bad"  points  for 
the  entire  interferogram,  it  was  ignored,  for  most  cases  the  interpolation  used 
added  only  a  minute  and  tolerable  changes  to  the  low  frequency  behavior  of  the 
interferogram.  In  the  case  that  one  of  the  extremums  (minimum  of  maximum)  of 
the  interferogram  is  found  to  be  a  "bad"  point,  tins  procedure  would  lead  to  an 
unacceptable  measurement  and  rejection  of  this  data.  However,  lor  all  the  cases 
checked,  the  extremums  of  the  interferogram  were  not  affected.  After  the  initial 
editing  was  performed,  the  program  logic  searched  for  the  center  of  the  inter¬ 
ferogram.  The  center  is  taken  as  the  maximum  excursion  of  the  interferogram. 

The  data  was  then  truncated  so  that  there  were  N  points  on  either  side  of  the 
center.  Thus  we  have  2N>1  points  in  a  final  edited  interferogram.  If  the  data 
did  not  have  N  points  on  either  side  of  the  center  of  the  interferogram  it  was 
eliminated.  This  algorithm  meant  that  we  could  compare  many  i n t e rferometer 
scans  and  be  assured  that  the  same  type  of  data  base  was  being  analyzed.  The 
resolution  frequency  of  the  data  is  improved  to  the  degree  the  truncated  inter- 
ferograms  occur  at  approximate  nulls.  To  the  degree  data  ends  at  nulls,  leakage 
(Fourier  Power  from  one  frequenev'  appears  at  neighboring  frequenc  i os )  is 
minimi  zed. 


1 ntor feromet or  Fourier  Analysis 

It  is  necessary  to  estimate  values  of  the  amplitude  and  phase  o  1*  the  luscrete 
Fourier  Transform  in  order  to  calculate  radiance  value-  from  the  i nt orforomet or 
mcasu remen ts . 


iis  t  imat :  ng  ampi  itude. 

i ho  resolution  frequency  of  tlu  results  of  a  last  Fourier  Transform  equals 

-,4*r  where  .J  is  the  number  of  points  used  in  the  last  lourier  Transform.  In  this 

*  T 

particular  Fast  lourier  Transform  the  number  of  points  used  was  J  =  2*4  (T  any 

integer).  To  dec  reuse  the  resolution  frequency  of  the  Fourier  transform  (not 

the  data)  and  to  use  the  fast  lourier  Transform  subroutine  as  efficiently  as 

possible,  the  edited  i nterferogram  measurements  were  continued  (padded)  with 

zeros  to  form  204 S  data  values.  The  actual  input  data  was  supplied  to  the  last 

Fourier  Transform  (I  FT)  as  follows. 


Let  RH  be  the  real  part  of  the  input  data  of  the  FFT  and  TIM  the  imaginary 

2 

part  of  the  input  data  which  is  assumed  complex  (Otnes  and  Iinochson*”)  .  Y.  is 
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'File  results  of  the  FIT,  after  unscrambling,  give  the  amplitude  and  phase  around 
t h e  m i  d  po  i n t  of  t  lie  inter fe ro g ram . 


;sjt  i mating  ])hase. 


From  the  results  of  the  FIT,  we  get  estimates  of  the  phases  of  an  inter¬ 
fere  gram.  If  the  Y  was  found  exactly  at  the  center  of  the  interferogram  the 
phase  estimates  should  be  either  0  or  >180  degrees.  Hie  measurements, 
however,  arc'  taken  at  a  spacing  of  /X  and  the  value  Y  occurred  somcwhcie  in 
tile  interval  -  -y  •  X  •  ^  where  X  is  the  offset  from  the  exact  center  of  the 
i n tor ferogram.  Thus  a  good  approximation  (under  the  assumption  that  the  trunca¬ 
tions  occurred  at  equal  distance  from  the  center  of  the  interferogram),  the 

.  ,,  .  .  ,  -  XT'*  1 ,80  W.* ]S0 

phase  will  hr  either  • -A-y  •  •  +  noise  or  -  •  -  +  180  +  noise,  where  .  is 

X  .-n 

the  t  requency  of  the  measurement  and  v.^  is  the  Nyquist  frequency.  ITic  analysis 
only  had  to  estimate  where  changes  in  the  above  phases  occurred.  Since  the 


Sb 


resolution  frequency  of  the  iuterferogram  was  quite  good,  the  phase  changes 
usually  occurred  between  consecutive  estimates  of  the  Hi.  Occasionally,  the 
changes  occurred  over  two  consecutive  estimates  of  the  phase. 

Background  Correction  ( instrument  funct  ion) 

The  amplitude  or  intensity  estimates  (A(u>)  l  of  an  actual  experiment  were 
corrected  for  the  intensity  inherent  in  the  instrument.  This  made  use  of  the 
knowledge  of  the  phase  estimates  (P(w))  from  the  actual  experiment.  The  procedure 
used  was  as  follows: 

1)  Find  B(oj)  which  are  estimates  of  the  amplitude  of  an  experiment  when  the 
onboard  interferometer  looked  at  a  cold  source. 

2)  The  values  of  P(m)  are  estimated  as  either  0  +  correction  degrees 
(call  0)  or  +180  +  correction  degrees  (call  +180).  Phase  P(u')  -  0 
degrees  was  assumed  for  frequencies  of  between  0  and  i\  .  P(ud  was  then 
assumed  to  be  +180  from  to  some  frequency  f2 .  It  was  tiien  assumed 
to  be  0  for  all  frequencies  above  f2 .  The  values  of  f \  and  fS  were 
found  by  searching  for  values  of  the  phase  where  |  PU+ZAi) -P(  *. )  j  ^  50 
degrees,  where  Aw  is  the  resolution  frequency  of  the  FIT.  When  P(i.) 
did  not  follow  the  expected  pattern  the  scan  was  ignored. 

3)  Form  corrected  estimates  C(m)  of  A(w)  either  C(u-.i  =  JA(.-)  -  B(w)  |  when 
P(m)  =  0  or  C(u>)  =  ( A(m)  +  B(w)  |  when  P(;o)  =  +180. 

Estimating  simulated  radiance. 

From  values  C(.o)  (corrected  intensity  values)  we  can  simulate  radiance 
values.  If  we  wish  to  simulate  the  discrete  values  of  a  radiometer  which  passes 
the  frequencies  between  ^  and  we  integrate  values  of  (.',(<;)  between  u-  and  j. 

Tli is  integral 


was  numerically  integrated  using  Simpson’s  Rule  and/or  a  trapazoidal  rule. 
Editing  simulated  radi anco  values. 

The  editing  of  the  simulated  radiance  used  a  similar  procedure  as  the  one 
used  to  edit  the  radiometer  data  with  two  changes.  The  first  change  wa s  that 
the  procedure  was  repeated  whenever  a  ’’bad”  point  was  detected  until  no  addi¬ 
tional  ’’bad”  points  were  detected.  The  second  change  was  that  whenever  an 


"event"  seemed  to  occur  after  the  initial  editing,  the  simulated  radiance  values 
were  processed  without  editing.  'Ill is  procedure  was  necessary  because  editing 
data  from  known  events,  although  it  loft  the  event  identifiable,  removed  some  of 
the  simulated  radiance  values  during  the  event. 

Con  cl  us i  on 

Interferometer  experiments  indicated  that  these  results  reproduce  the 
measurements  which  would  be  obtained  by  a  radiometer  with  any  particular  passhand. 
This  interferometer  technique  has  the  enormous  advantage  that  any  band  may  be 
selected  via  computer  processing  to  provide  information  equivalent  to  several 
hundred  individual  radiometer  measurements.  The  timing  and  positioning  of  the 
experiment  described,  however,  was  different  according  to  which  direction  the 
interfere gram  was  sampled.  This  led  to  Fourier  Power  near  the  Nyquist  Frequency 
of  the  computed  radiance.  If  this  Fourier  Power  was  ignored  the  resulting  radiance 
values  led  to  results  which  were  consistent  with  the  results  of  radiometer  experi¬ 
ments.  Other  problems  with  this  technique  are  that  the  analysis  is  not  straight¬ 
forward  and  a  very  large  data  base  is  needed  for  the  computations. 

epilogue 

Further  relevant  discussion  and  data  analyses  procedures  will  be  presented 
in  an  AF(.IL  technical  memorandum  dealing  with  a  BAMM  5  data  report. 
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PRIMARY  AURORAL  ELECTRON  PRECIPITATION  PROBLEM  AND 
ELECTRON  AND  PARTICLE  DISTRIBUTION  FUNCTIONS 

Over  the  last  two  years  mathematical  and  numerical  analysis  techniques 

have  been  applied  to  describe  the  effects  of  electrons  in  the  ionosphere.  In 

particular,  two  types  of  problems  were  studied.  The  electron  and  particle 

distribution  functions  in  the  ionosphere  were  predicted  using  computer  code 

described  In  the  Report  AFGL-TR- 79-01 32 .  During  the  last  two  years  this  code 

was  modified  to  incorporate  the  latest  models  for  the  collision  and  energy 

transfer  cross  sections  for  the  particles  in  the  ionosphere.  In  particular, 

referencing  the  report  AFGL-TR- 79-0132,  Q1 .  (E.)  is  now  calculated  using  either 

J  >  ^  * 

of  two  formulas.  The  first  formula  states  that 


if  El.  .  >  E.  orWL.  _ 
3,l~i  J,L 


li; ,  then  QI  „([;.)  =  0 


Otherwise 


81  , 
AI.  v  J,L 
J.l*  \ 


where 


WI.  y 

J,L 

E, 

1  - 


01.  ,  if  01.  . 
3,1-  3,L 


C  .(°h±h\ .  2 

V  \  WIJ.L  ) 


718281828 


if  01 


Alternately  the  second  formula  states  that  QI 


i.ii'V 


0  if  11.  , 
J  ,L 


Otherwise  QIj-L(E.J  -  1016- (KI  j  L/ Uf )  wbN^y-i-^*  I  j  L*li.  /  C  U.  .G«  1 .  ( J ) 
(tan1 ( (. 5*(E. -FI  ) - (TS I  .-(TAT  /(E . *TBI  .))))/ 

1  J  i  o  J  *  •  •  .1  »  l J  •  J  t 
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( t  ( IS  I  *t.  J/(n.+C»lj(L)))*tan_1(TSI  iL-(TAl  y  ll-.+TBI  .  L 

*  (  I:,  +(ilil  .  ,*h.)) 

1  J  y  Ij  J  y  \a  J 


J)J 


In  these  representations  El  .  ,  ,  CS I. 


KI 


hi  i  iiv.tv  i>  i  .  .  ,  ixi  i  •  ,  t  QO  ,  GB I .  .  ,  TB I .  .  ,  TS 1  .  .  , 

J  >  I'  J  >  L,  J  >  *-■  J  >  U 

,  ,  JI.  ,  ,  TAI .  .  ,  HII.  .  ,  El.  .  ,  01.  .  ,  AI.  .  ,and  BI.  .  are  input  data 

L  J,L  ij  i  ,;i  i,j  i,j  i,J  J,L 

IO  Q 


val ues 


Another  modi f icati on  was  a  redefin  ition  of  Q„  ,,  in  particular  for  T=l,  2 

l ,  K 

or  3.  'lhc  new  formula  takes  the  form 


V  “°  if  Ei  i°T.K 


Otherwise  Q„, 

i  ,  K 


<>0^|-.k/<11J.K2»*FX' 


where 


LX  = 


U 


T,  K 


E. 

1 


T,  K 


1 f  °T,K  '  ‘  7 


and 


LX  = 


HT  .A  /  0  *E.  \ 

*  LN  4.  *  — — -  +  2.718281828 

h  )  \  "t,k  / 


if  °T,K  —  '7 

In  addition  some  values  of  (j  which  were  originally  found  by  formulas 

I  ,  K 

are  now  supplied  as  tables  of  values  versus  IE.  These  and  other  changes  were 
made  to  the  computer  model.  Drum  plots  were  also  made  for  all  the  cross-sections 
now  used  by  this  model. 

Hie  second  major  effort  performed  for  this  problem  was  calculating  electron 
and/or  proton  fluxes  in  auroral  el ectron/proton  precipitation.  Several  mathe¬ 
matical  models  describing  this  problem  have  been  studied.  The  following  discus¬ 
sion  describes  the  model  upon  which  greatest  efforts  were  expended  to  determine 
the  el  eel  ron  flux. 


Let  this  electron  flux  be  <J>(t,L,vO  where  t  is  optical  depth;  E  is  the  energy 
and  \i  is  the  cosine  of  the  zenith  angle.  Then  a  trans{>ort  equation  describing 
the  electrons  is 
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♦  O.E.u)  =  *(E)  eT/tJ  -  J 


d_t  e  (T-t)/u 
M 


b(E)  4>(b(n)  *t,E+W,yJ  -  r(E)*<|»(t,H,y) 
r(E)  J  dy'  p(E,  y* ,  y)  *4>  Ct  ,H,  y’)j 


-1  <  n  <0  and 


o 

<Kt,E,u)  =  J 


dt  ^(T-tJ/u 
V 


jb(E)*«j>(b(E)*t,  E  +W ,  y)  -  r(E)*4>(t,n,p) 

rl  1 

r(E)  /  dy '  p(E,y',y)  y'  )> 

-1 

d 


for  0<  p < 1 . 


/KJQs  \ 

■(jrwEcxp  t_  c/,y 

and  K, ,  Q  ,  and  E  are  constants  with 
1  o 

b CE)  =  21^1  ,  r(E)  -  QE'  (K)/Q(ll) 

Q(EJ  and  QE’ (E)  arc  defined  by  logarithmic  interpolation  o 
versus  E  and  QE* (E)  versus  E  respectively. 


1)  1  os  o f  ()(  l: ) 


P(E, , „) . _ 1 1 ; 

f(l»:n(E)-u'u)2  -  (l-u"Vll-P-Jl] 

where  n(E)  is  found  by  logarithmic  interpo  1  at  ion  of  a  table  of  n(E)  versus  I: 


values , 


The  integral  equat ion s  were  approximated  for  a  set  ot  N^E.  values  I*  ^  , 

tMTkI  +  W,  Emtm  +  2W,  .  .  .  , H...  a  set  of  N  ,t.  values  0,  i  .  ...  t  (l  • 

MTN  MAX  i  j  2  MAX  v 


and  a  set  of  N  values 
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T he  p  integrations  ^ J  dp*  etc^  wlto  rewritten  as  ^ ^  dp'  etc^  + 


^  J  dp1  etc.^  bcc:j 


:ausc  nici y  be  discontinuous  at  <f>(r,E,  o) .  The  following 


procedure  describes  t  ho  technique  used  to  determine  <j>  (*u ,  i.,o )  . 


Rewriting  ^  1 21  and  (.13),  we  have 


h,  v  ^  -if\\  /  dt  ( T  “  t )  /  p 

-  4^ C *  )  O  -  /  —  e 

Jo  U 


for  i  •  u  <  0  and 


>H  t  >  *•»  M) 


=  fd±  c(l't,/qru,i:.1i)) 


where 


l?U  .Ii.  m)  =  b(l-)  *^(b(E)*t,l:+VV,ii)  -  (r(E)  *•;  (t , E,  p)  +  r(E) 


J  dp*  p(i:,p\p)*<|>(t,E,p')  . 


As  n*>0“  in  (16)  and  p  >0+  in  (17),  F(t,E,p)=>  F(t,I‘,0”)  and  F(t,E,0+) 


respect i vel y . 


Thus  cf> ( t 


t 

.  b,  ()*)  =  -!••(!, K, O')  J 


dt  (T-t)/0‘ 


etT'tj/0’  *  0  for  to’ 


Hence, 


)  -FOJ 


S i mi  1 ar 1 y 


:,0')  [ 


.(  t-t)/0" 


Ffi,r.,o')  . 


{'(  r ,  E,0+J  -  +l;(  i ,  l.,0+  )  /  ^(  t  -tj/0  =+f(t>,:>0  +  )  J 


Tr(T-t)/cr 


-  K  t.r.,o  )  . 

Now  consider 

<;(  i,i:,m)  --  i(i  )  N'd  a.,i  )  +  r(\-)*  J  dp'  p(i:,ii',p)  *<j.(t,E,pS 


with  the  assumption  that  for  values  of  p'  such  that  p  <  u*  <  u 

k  k  + 1 


QCt.E.p*)  = 


tlJk+rp,) 


a'k+ruk} 


)  (  W  1  ~  M,  ) 

j  *  r„k— 


where  linear  Interpolation  is  used.  Ihcn 

v*  /r" 

=  -r(li)  +  r(i:j  £  (  I 

k=l  V^u, 


ill''  p(i;,l»\u)*<Kt,l:,u,J 


where 


du'  p(L',u',y)  =  0 


vyfv, 

^•G(t,E,u)  =  -r(E)  <Kt,p,E)  +  r  (K)  Eli  <Jp'  pdi.u'u) 

k*l\X 


Mk+1 

- - -  <t>(t,E,p  )  -  - - —  )  -  - - --  <t> C t .  l;» 

Mk  +  1  \  k  "k  +  l  *  Mk  k  ,:k  +  l  '  lJk 


-  u.  ^  ’  ’  k+r 


k  +  l  Mk 


Now  consider 


p(E,p',p)  = 


2n(E)  *(  l+n(E) )  *  ( 1  >i  ( ti)  - 

”  ~~  *7>/  f 

J(l+2ri(HJ -p(j'  )*■  -  ( 1  -  n  •  *" )  *(.  l-p")J 


a  2n(Ej*(1 +n(E)  )*(I+2ri(E))  -  2n(F.m*n(U)  *m>’ _ 

pri(E)  +  4n“(E)  +  u“  -  (2p+4n(Ej  *i.)  |i *  ■<  (:’“J 

_ _ u_ +_\v _ 

r  2 

p+bp'+ii' 


in  <MjuntioJi  (24), 


u  =  L’n(io*(i+n(E))*(i*2n(i;)j,  v  =  -2n(i:)*(i+n(t))*iJ  ) 

>  b  =  -(2y+4n(H)*„) 


I 


a  -  4 n ( 1  )  +  Iff  (l:j  +  u 
'/here fore,  if  u  f  +I 


K+1  J  M  f 

du’  p(n.M’M)  (-"■■_  ■  Ht.b.v .)  -  --- II...  »(t>E,p  ) 

Mk  +  1  pk  K  k 


"k  +  1 


uk+l  llk 


u' 


k+1'  "  Cm^.-iv)  ♦(t»E»Mk+l) 


r  4 

-'ll.  \ 


k+1  *k 

A,i  ♦  Vt »’  *  ck,k"' 


>) 


3/2 

R  '  (u') 


<Kt,E,pk) 


'k,k+i  Bk,k+i  u  Lk,k+1  M 

3/2 

R  O’) 


Kt,E,uk+1) 


where  Kip’)  =  ri  +  bp'  +  p'' 


k+1 


k,k  ‘-vrV  ’  k,k 


■  v  "  tv  uk  +  1-.))/(pk+1-pk)  , 


ck,k  =  'v/("k+r|ik)’  \ , k + 1  =  *u  Tuk+1  -pk)  •  and 


*k.k+i  =  (u‘Vpkl/(“k+rv,k)> 


-k.k+j  o.k+rukj  • 


f.ct  t  i n  u 


a  --  ia-b2  -  lOn(H)*(j*n(U)rU-n2)  ,  A/o 


Now  p  will  be  evaluated  at  one  of  t  ne  set  of  p_,  values  fsav  p  )  and 

m  k  v  -  nr 

i  N 

f  ,J 

-rOO  +  rcii)  J  d;:'  P(l.  f* ,  =  £  *1Kflini)  (3Jj 


.)  -  r (. Li)  CD,  .  +  1),  ,  .  -  rt  .) 

k  k ,  k  k  - 1  ,  k  m ,  k 


w,rh  %.l  =  “n  ,n  r  0 

M  P 


I  f  ;•  -  p  ;r  then 

in 


P^pSpO  = 


--~2 
ld+ep'  ) 


where  C=2n(h)  U  +  n(L) ) ;  d^l+Jn(h)  and  e  -  -  p  •  liquation  (2(>)  for  p  =  p  =  +1  now  becomes 

m  m 

rM|<+l  /  mk+i  » 

du'  p(li,  l*  *  >  pi  (f~  r  <p(t,i:,p)  -  7—--7T  4>(t,E,p  ) 

Jpk  Y'  K+l  "kJ  K  (MK+rMk}  K 


fv?"1'''"1"1  *  1^5?  ♦(,-i"W 


i:.  ,  l:,  ,  p' 

♦it.n.iO  +  — - — <Kt.E 

(.il+t  K  tU+o 


(  k  ,  k  + 1  ,  .  k , k  +  l 

Itt.t-  *‘,,Ii’,ikm)  f - 

\1  d  +r  *  (i '  )  (dU-*,i')‘ 


w  he  re 


r;  .  r '  fk  >1  _  -1; 

k,k  (UK+1-UK)  ’  K.K  ”  (mk+1-VkT 


k  ,k+1  ^K+V^d  KK+‘  i;k+rT;K 


However, 


f  du' 

J  (d+e  *u 


e(  d+e*|j '  J 


f  du1  u* 

jcd+cTT2 


— — - +  iLn  (,d+e*u'j 

e(d+e*M’)  e 


Therefore  for  u  =  u  =  +1 , 

ni 


HK,K  +  t;K,K 

n 

( d+e *11 '  J”  ' 


4>(t ,  1-,  Uj-3 


/^LK,K+1  1  K.K-.  1^\  x,x 

+  (  — 2 - ^ - j  J 

\  (d+c*M  * } 2  /  K  ‘ 


=  Dk)k  Ht,K,VK)  ♦  uk)k  +  1  <Ht.U,PKtl) 


When  =  pt  =  +1 , 


D  -  I  A_i _ 1 _ \  +  [  I  f  _ A+L_ 

k,L  e  '^(d+e*MK)  (d+e*PK+1W  K,L  l  e  ^(d+oAMK)  (<l+c*’!  K  +  ly 


i  /d+c*Vi\\, 

+  e2  *"(d+CX  JJ'-K.I.  * 

For  L=k  and  5,  =  k+1  in  the  above  expressions,  (SI)  then  holds  lor  .  * 

Now  consider  the  term 
T  =  bCE)**(b(E)t>E+W>gm) 

We  wish  to  evaluate  T  for  !:=!£.,  i~l,...,N^.  I  f  ^  \v|\x  ‘  »  1  ';,n  lH* 

linearly  interpolating  on  t  to  £Ct 


found  hv 


T  =  b(E) 


/O^j-bUOt) 

0 


(b  t  in t - 1  ) 

(T  .  -1  .  ) 

J  + 1  I 


'Hr.  .  Ji+Wpi  j1 
i  +  1  in ; 


where  and  H<h1(.v  .  Then 

.1-  —  j+i  MAX 


T  =  II.  .  *<{>(t  •  ,l-fW,p  J  +  1  •  *t  *^(x  .  ,i:+K,  ji  ) 

J .  .1  J  m  1,1  .1  inv 

+  u.  .  ’-ht.  -,n+w,ii  )  +  iL  .  *t * ; ( ; .  ,  ,l:+w, i>  ) 

J.J+1  j  +  1  nr  j,j  +  1  j+1  nr 

who  I'  l' 


j  MH)  t.  , 

r  -  L+J 

i .  i  It.  ,  -x  . ) 
J+1  J 


h(h)  7 

IV  =.  - - J.. 

J  ,  I  +  1  t  T  .  .  -  l  .  ) 

J+1  I 


r:  .  = 


_!» 2_m_ 


I  ,  J  T  .  .  -  T  . 

J+1  | 


,L  -  b-(E) 

'•J*1  \,.i  -  ’j 


If  h.  *  ^M/W  =  lij\'  }  t  {*  assumed  that  as  !.-**> 


-ot(x . , u  )  W 

*Y’L+,V’V  -  '  m 


where  a(t.,u  )  is  ,i  set  of  input  data  values.  lor  the  ease  H=EW, „  we  will  then 
J  111  MAX 


T  =  "jJX  ‘*(YW>V>  +  ‘ij"  *^YWV 


J  .  ,i  + 1  '.i+ 1  ,l'MAX,;ini)  *  ’  j  ,  j+1  lj  +  r,:MAX,w»J 


ii  MAX  _  *  ,(I'MAX)  li  ‘'MAX  „  *  b(l:MAX^  Th 

"...I.  " k  a,,J  ha.  ■ k  w  h,i. 

In  Inc  above  equation,  I;  1 


f « 1 1*  i.  j  and  i.  ~  j  *  1  ,  and  K 


'max’*  al  VV  w 


Isinj4  the  above  result  *  we  now  wish  to  calculate 


IK 


S^t.K.u)  =  /  ^  e(T_t)/u 


{l'(t 


s  tCt.e,  u) 


rw 

K  v 


i-t)/u 


From  equations  (19)  and  (20) 
=  F ( t , H , o )  =  S2 


When  u  ?  o  ,  u  =  p  , 
in 

EW  =  E+W*(1-6C  r  ) 


h(Fj  is  assumed  <1,  •*'  ,  and  x  ^  <  t  •-  r  .  These  assutnpt  ions  r;m  J>e 

satisfied  so  long  as  !<!,.. 

&  —  Nt 

We  can  then  express 
F(t>E’pid  =  Hj,j  *^Tj’EW>^  + 

. ~'^{L*l**(Tj  +  l,CW’,Jin)  +  >  J  +  1  *t*^(T  j  + 1  ’  l'*v ’ "m^ 


k J  lJlyl 


K«.,£+l  *^Tfc+l,E,,Jld  +  lSL,  M*1  **^Ti+l'L,,lkJ  ' 


where  wc  have  used  equations  (311,  (*U)  and  (431  while  assuming  linear  i  nt  e  rpn  1  at  i  on 


between  xf  and  ..  Now 

l  P.+ 1 


'  *•»  (Vi-v 


K  -  1  I  - _ L _ 

*-t+1  "  ’  <WV  ’  '^+r  (TurT 


I'ii  order  to  evaluate  in  ( 44 J ,  we  note  that 


C  'i  •  i  •  ,) 
.1  J  - 1 


Sjltj  ,l:,n)  =  C 


dt  ^(x  tj/u 

U 


.  (4i>J 


ing  tiiis  fact  wo  first  evaluate  S^(t9,E,p)  =  0  since  t^=0) 


We 


then  use  this  result  to  evaluate  S^(i_.,E,pJ  and  continue  this  procedure  recursive¬ 
ly  until  we  evaluate  S  (t^  ,  E,p). 

T 

Similarly  in  order  to  evaluate  S7  in  (44)  ,  we  note  that 

SUij.L, MJ  -  e(  J  J+l)/'  SJ(Tj  +  1,li,p )  *J  j  ~  e(T't)/M{F(t,E,p)}  .  (50A.) 

1  j 

Tii us,  we  first  evaluate  S^(tN<  ,  J:,p).  We  then  use  this  result  to  evaluate 

_  N 

T 

S^(t  and  continue  this  procedure  recursively  until  we  evaluate 

^(^=0,1:, PJ. 

Enough  information  has  been  provided  to  calculate  the  and  S0  values,  so 
long  as  S((t^  ,li,p)  can  i>e  evaluated.  In  order  to  evaluate  S ,E,p),  the 

I  T 

fo  1  1  ow  j  n i\  mode  1  i  s  used  : 


■*et  pn  =  a  and  i»(N  )  ^  p 


(SOB  J 


-  a(E,  n  .  ) 


iMii  Lh  1. .  .  llj_ 

1  M  —  I  1  ~ 


The  values  of  are  found  by  interpolating  a  table  of  t(E,u)  versus  I.  and 
M,  which  is  provided  as  input  data.  The  values  of  $(T,litu)  with  r>;^  are  then 
assumed  to  be 


-  (  i  i  v  )  ot(  E,  \i) 

i\ 


C'i  ,  m  uj  -  <-• 


SO 


With  this  model  we  calculate  S0(a^  ,li,P) 

T 


f~  <’n  -t)/u» 

S2tCN  -J  rC  ’ 

T  tn  m 

T 


-  -  f ' 


dt  (Vt)/w« 

—  e 


bxXT  m 
N 

T 


hi*T  .  <Kt  -,EW,U  )  +  1*.’  .  *t  *  ,'■( .  .  ,  El\, ;.  ) 
(_  J ,  J  J  nr  j,.l  '  111 


*6(t  EW  u  1  +  T  *t*'hfT.  .  , l.W ,  i i  ) 

Hj,j+1  j+1*  'V  j ,  J  +  l  V  ,l+l 


^nt  ~  h  (E) ' 
ds 


Cs-tn  )  a(E.,Um) 


♦lTN  'EN»wm)  ° 

T 


-  f  dt  <G(t ,  E,  u^)  e 

C 


where  the  integration  from  bt^  to  tn  is  described  in  (531  and  (54).  The 

T  T 

integral  over  ds  is  found  by  the  substitution  s=b(EJt  and  is  described  in  (55) 
and  (56).  The  integral  of  C.(c, H,n  )  (see  equation  (21))  appears  in  (64). 

Thus  we  have 


(t-t) 

il  e  Um 


n 


£'L  *  'j,j 


<V>  ,f  (T-t )/,  f  ,  ,  . 

V  /  c  hiL  .  <j)(T . ,  ew.  i.  )  +  i  v '  .  t 

L  /  u  /  i »  i  »  m 

j=l  -'as  )  %  L 


("Li 


r  m 


, EW .a  )  1 
m  ; 


who  iv 
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Using  the  preceding  npproximnt  ions  for  the  behavior  of  <J>(t,E,u)  and  the 
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This  compl ot  cs  tlu>  doscript  ion  of  methods  used  to  determine  the  electron 

flux 


